function [modelSim,moments] = simulateModel(xGrid,xDrift,xVolatility,modelSolution,algParams,pruneXSim)
%SIMULATEMODEL Simulates the model
%   Inputs:
%     - xGrid: grid for state variable x (in practice, x=\tilde{\sigma}^2,
%              but this is not assumed)
%     - xDrift: (arithmetic) drift of state variable x on the grid
%     - xVolatility: (arithmetic) volatility of state variable x on the grid
%     - modelSolution: struct array containing the model solution (output
%                      of function solveModelPde)
%     - algParams: struct array containing algorithm parameters (requires
%                  the parameters nSimPeriods, nSimBurnIn, simTimeStep)
%     - pruneXSim: flag indicating whether x dynamic should be pruned to
%                  ensure that x never leaves [min(xGrid),max(xGrid)]
%   Outputs:
%     - modelSim: struct array containing simulated paths
%     - moments: struct array containing computed model moments


  % unpack modelSolution
  [sigmaI,a,govSpending,theta,growthRate,consCapRatio,qB,qK,iota,equityPremium,equityExcRetVol,adjBondGrowth,growthRateCons,riskFreeRate,expectedReturnBonds,sigmaTotalWealth] = ...
    aux.unpack(modelSolution,'sigmaI','a','govSpending','theta','growthRate','consCapRatio','qB','qK','iota','equityPremium','equityExcRetVol','adjBondGrowth','growthRateCons','riskFreeRate','expectedReturnBonds','sigmaTotalWealth');

  % unpack algParams
  [nSimPeriods,nSimBurnIn,simTimeStep] = ...
    aux.unpack(algParams,'nSimPeriods','nSimBurnIn','simTimeStep');

  stepsPerPeriod = 1/simTimeStep;
  simShocksBurnIn = randn(nSimBurnIn*stepsPerPeriod,1);
  simShocks = randn(nSimPeriods*stepsPerPeriod,1);

  stateDriftInterpolant = griddedInterpolant(xGrid,xDrift);
  stateVolatilityInterpolant = griddedInterpolant(xGrid,xVolatility);
  growthRateInterpolant = griddedInterpolant(xGrid,growthRate);
  
  if pruneXSim
    xMin = min(xGrid);
    xMax = max(xGrid);
  else
    xMin = -Inf;
    xMax = Inf;
  end

  [~,minIdx] = min(abs(xDrift));%find stochastic steady state
  xStart = xGrid(minIdx);
  xBurnInState = xStart;
  for i=1:(nSimBurnIn*stepsPerPeriod)
    xBurnInState = xBurnInState + stateDriftInterpolant(xBurnInState)*simTimeStep ...
      + stateVolatilityInterpolant(xBurnInState).*sqrt(simTimeStep)*simShocksBurnIn(i);
    xBurnInState = min(max(xBurnInState,xMin),xMax); %prune simulated state
  end

  simX = NaN(nSimPeriods*stepsPerPeriod+1,1);
  simX(1) = xBurnInState;
  simK = NaN(nSimPeriods*stepsPerPeriod+1,1);
  simK(1) = 1;
  simTime = NaN(nSimPeriods*stepsPerPeriod+1,1);
  simTime(1) = 0;
  for i=1:(nSimPeriods*stepsPerPeriod)
    simX(i+1) = simX(i) + stateDriftInterpolant(simX(i))*simTimeStep ...
      + stateVolatilityInterpolant(simX(i)).*sqrt(simTimeStep)*simShocks(i);
    simX(i+1) = min(max(simX(i+1),xMin),xMax); %prune simulated state
    simK(i+1) = simK(i) + simK(i)*growthRateInterpolant(simX(i))*simTimeStep;
    simTime(i+1) = simTime(i) + simTimeStep;
  end

  % compute some relevant model variables
  sigmaInterpolant = griddedInterpolant(xGrid,sigmaI);
  equityPremiumInterpolant = griddedInterpolant(xGrid,equityPremium);
  equityExcRetVolInterpolant = griddedInterpolant(xGrid,equityExcRetVol);
  thetaInterpolant = griddedInterpolant(xGrid,theta);
  consCapRatioInterpolant = griddedInterpolant(xGrid,consCapRatio);
  aInterpolant = griddedInterpolant(xGrid,a);
  qBInterpolant = griddedInterpolant(xGrid,qB);
  qKInterpolant = griddedInterpolant(xGrid,qK);
  iotaInterpolant = griddedInterpolant(xGrid,iota);
  adjBondGrowthInterpolant = griddedInterpolant(xGrid,adjBondGrowth);
  govSpendingInterpolant = griddedInterpolant(xGrid,govSpending);
  growthRateConsInterpolant = griddedInterpolant(xGrid,growthRateCons);
  riskFreeRateInterpolant = griddedInterpolant(xGrid,riskFreeRate);
  expectedReturnBondsInterpolant = griddedInterpolant(xGrid,expectedReturnBonds);
  sigmaTotalWealthInterpolant = griddedInterpolant(xGrid,sigmaTotalWealth);

  simTheta = thetaInterpolant(simX);
  simSigmaI = sigmaInterpolant(simX);
  simEquityPremium = equityPremiumInterpolant(simX);
  simEquityExcRetVol = equityExcRetVolInterpolant(simX);
  simConsCapRatio = consCapRatioInterpolant(simX);
  simConsumption = simConsCapRatio.*simK;
  simA = aInterpolant(simX);
  simOutput = simA.*simK;
  simQB = qBInterpolant(simX);
  simQK = qKInterpolant(simX);
  simIota = iotaInterpolant(simX);
  simInvestment = simIota.*simK;
  simGovDebt = simQB.*simK;
  simAdjBondGrowth = adjBondGrowthInterpolant(simX);
  simPrimarySurplus = - simAdjBondGrowth.*simQB.*simK;
  simGovSpending = govSpendingInterpolant(simX).*simK;
  simGrowthRateCons = growthRateConsInterpolant(simX);
  simRiskFreeRate = riskFreeRateInterpolant(simX);
  simExpectedReturnBonds = expectedReturnBondsInterpolant(simX);
  simSigmaTotalWealth = sigmaTotalWealthInterpolant(simX);

  % select quarterly data
  offset = 0.25/simTimeStep;
  simOutputQuarterly = simOutput(1:offset:end);
  simConsumptionQuarterly = simConsumption(1:offset:end);
  simSurplusQuarterly = simPrimarySurplus(1:offset:end);
  simSurplusOutputRatioQuarterly = simSurplusQuarterly./simOutputQuarterly;
  simInvestmentQuarterly = simInvestment(1:offset:end);
  simGovSpendingQuarterly = simGovSpending(1:offset:end);
  simGovDebtQuarterly = simGovDebt(1:offset:end);
  simGovDebtOutputRatioQuarterly = simGovDebtQuarterly./simOutputQuarterly;

  [~, outputLogHP] = hpfilter(log(simOutputQuarterly),1600);
  [~, consLogHP] = hpfilter(log(simConsumptionQuarterly),1600);
  [~, surplusOutputRatioHP] = hpfilter(simSurplusOutputRatioQuarterly,1600);
  [~, invLogHP] = hpfilter(log(simInvestmentQuarterly),1600);
  [~, govSpendingLogHP] = hpfilter(log(eps+simGovSpendingQuarterly),1600); %add eps to accommodate a solution without gov spending
  [~, govDebtOutputRatioHP] = hpfilter(simGovDebtOutputRatioQuarterly,1600);

  % compute moments
  stdOutput = std(outputLogHP);
  stdCons = std(consLogHP);
  stdInv = std(invLogHP);
  stdGovSpending = std(govSpendingLogHP);
  stdSurplusOutputRatio = std(surplusOutputRatioHP);
  acorrOutput = corr(outputLogHP(1:end-1),outputLogHP(2:end));
  corrSurplusOutputRatioOutput = corr(outputLogHP,surplusOutputRatioHP);
  corrConsOutput = corr(outputLogHP,consLogHP);
  corrInvOutput = corr(outputLogHP,invLogHP);
  corrGovSpendingOutput = corr(outputLogHP,consLogHP);
  meanConsOutpRatio = mean(simConsumption./simOutput);
  meanInvOutpRatio = mean(simInvestment./simOutput);
  meanGovSpendingOutpRatio = mean(simGovSpending./simOutput);
  meanTheta = mean(simTheta);
  meanCapitalOutpRatio = mean(simQK.*simK./simOutput);
  meanDebtOutpRatio = mean(simGovDebt./simOutput);
  meanSurplusOutputRatio = mean(simPrimarySurplus./simOutput);
  meanEquityPremium = mean(simEquityPremium);
  sharpeRatio = mean(simEquityPremium)/mean(abs(simEquityExcRetVol));
  meanSigmaI = mean(simSigmaI);
  stdSigmaI = std(simSigmaI);
  meanIota = mean(simIota);
  meanAdjBondGrowth = mean(simAdjBondGrowth);
  meanRiskFreeRate = mean(simRiskFreeRate);
  stdRiskFreeRate = std(simRiskFreeRate);
  stdDebtOutputRatio = std(govDebtOutputRatioHP);

  % create outputs
  modelSim = createSimOutputStruct(simTime,simX,simSigmaI,simA,simIota,simConsCapRatio,simTheta,simQB,simQK,simEquityPremium,simEquityExcRetVol,simAdjBondGrowth,simK,simOutput,simConsumption,simGovSpending,simPrimarySurplus,simGovDebt,simGrowthRateCons,simRiskFreeRate,simExpectedReturnBonds,simSigmaTotalWealth);
  moments = aux.pack(stdOutput,stdCons,stdInv,stdGovSpending,stdSurplusOutputRatio,acorrOutput,corrSurplusOutputRatioOutput,corrConsOutput,corrInvOutput,corrGovSpendingOutput,meanConsOutpRatio,meanInvOutpRatio,meanGovSpendingOutpRatio,meanTheta,meanCapitalOutpRatio,meanDebtOutpRatio,meanSurplusOutputRatio,meanEquityPremium,sharpeRatio,meanSigmaI,stdSigmaI,meanIota,meanAdjBondGrowth,meanRiskFreeRate,stdRiskFreeRate,stdDebtOutputRatio);
end

function modelSim = createSimOutputStruct(time,x,sigmaI,a,iota,consCapRatio,theta,qB,qK,equityPremium,equityExcRetVol,adjBondGrowth,k,output,consumption,govSpending,primSurplus,govDebt,growthRateCons,riskFreeRate,expectedReturnBonds,sigmaTotalWealth)
  modelSim = aux.pack(time,x,sigmaI,a,iota,consCapRatio,theta,qB,qK,equityPremium,equityExcRetVol,adjBondGrowth,k,output,consumption,govSpending,primSurplus,govDebt,growthRateCons,riskFreeRate,expectedReturnBonds,sigmaTotalWealth);
end
